Overlays

HES 505 Fall 2022: Session 18

Matt Williamson

Objectives

By the end of today you should be able to:

  • integrate a covariate into KDE’s

  • Describe the utility and shortcomings of overlay analysis

  • Describe and implement different overlay approaches

Re-visiting Density Methods

  • The overall intensity of a point pattern is a crude density estimate

\[ \begin{equation} \hat{\lambda} = \frac{\#(S \in A )}{a} \end{equation} \]

  • Local density = quadrat counts
nclust <- function(x0, y0, radius, n) {
              return(runifdisc(n, radius, centre=c(x0, y0)))
            }
x <- rpoispp(lambda =50)
Q = quadratcount(x)
plot(x, main="")
plot(Q, add=TRUE)

Kernel Density Estimates (KDE)

\[ \begin{equation} \hat{f}(x) = \frac{1}{nh_xh_y} \sum_{i=1}^n k\bigg(\frac{{x-x_i}}{h_x},\frac{{y-y_i}}{h_y} \bigg) \end{equation} \]

  • Assume each location in \(\mathbf{s_i}\) drawn from unknown distribution

  • Distribution has probability density \(f(\mathbf{x})\)

  • Estimate \(f(\mathbf{x})\) by averaging probability “bumps” around each location

  • Need different object types for most operations in R (as.ppp)

Kernel Density Estimates (KDE)

  • \(h\) is the bandwidth and \(k\) is the kernel

  • We can use stats::density to explore

  • kernel: defines the shape, size, and weight assigned to observations in the window

  • bandwidth often assigned based on distance from the window center

x <- rpoispp(lambda =500)
K1 <- density(x, bw=2)
K2 <- density(x, bw=100)
K3 <- density(x, bw=2, kernel="disc")
Warning in density.ppp(x, bw = 2, kernel = "disc"): Bandwidth selection will be
based on Gaussian kernel
par(mfrow=c(1,3))
plot(K1, main="Bandwidth=2")
plot(K2, main="Bandwidth=100")
plot(K3, main="Bandwidth=2, kernel=disc")
par(mfrow=c(1,1))

Kernel Density Estimates (KDE)

  • What if we imagine intensity is a function of another covariate
load(url("https://github.com/mgimond/Spatial/raw/main/Data/ppa.RData"))
ma <- readRDS("file/ma.rds")


# Load an MA.shp polygon shapefile 
    s  <- as(ma, "sf")
    w  <- as.owin(s)
    w.km <- rescale(w, 1000) 

    # Load a starbucks.shp point feature shapefile
    s  <- st_read("file/Starbucks.shp", quiet=TRUE)  
    starbucks  <- as.ppp(s)
Warning in as.ppp.sf(s): only first attribute column is used for marks
    marks(starbucks) <- NULL
    Window(starbucks) <- w     
    
    # Load a pop_sqmile.img population density raster layer
    img  <- raster::raster("file/pop_sqmile.img")
    pop  <- as.im(img)
Q <- quadratcount(starbucks, nx= 6, ny=3)
plot(starbucks, pch=20, cols="grey70", main=NULL)  # Plot points
plot(Q, add=TRUE)  # Add quadrat grid    

plot(log(pop))

Adding a Covariate to KDE

  • rhohat computes nonparametric intensity \(\rho\) as a function of a covariate

\[ \begin{equation} \lambda(u) = \rho(Z(u)) \end{equation} \]

starbucks.km <- rescale(starbucks, 1000, "km")
ma.km <- rescale(w, 1000, "km")
pop.km    <- rescale(pop, 1000, "km")
pop.lg <- log(pop)
pop.lg.km <- rescale(pop.lg, 1000, "km")
rho <- rhohat(starbucks.km, pop.lg.km,  method="ratio")

Adding a Covariate to KDE

plot(rho, las=1, main=NULL, legendargs=list(cex=0.8, xpd=TRUE, inset=c(1.01, 0) ))
Warning in min(D[scaledlegbox]): no non-missing arguments to min; returning Inf

pred <- predict(rho)
cl   <- interp.colours(c("lightyellow", "orange" ,"red"), 100) # Create color scheme
plot(pred, col=cl, las=1, main=NULL, gamma = 0.25)

Adding a Covariate to KDE

  • We can also think more generatively

  • Model explicitly as a Poisson Point Process using ppm

\[ \begin{equation} \lambda(u) = \exp^{Int + \beta X} \end{equation} \]

# Create the Poisson point process model
PPM1 <- ppm(starbucks.km ~ pop.lg.km)
Warning: Values of the covariate 'pop.lg.km' were NA or undefined at 0.57% (4
out of 699) of the quadrature points. Occurred while executing: ppm.ppp(Q =
starbucks.km, trend = ~pop.lg.km, data = NULL, interaction = NULL)
# Plot the relationship

Adding a Covariate to KDE

PPM1
Nonstationary Poisson process

Log intensity:  ~pop.lg.km

Fitted trend coefficients:
(Intercept)   pop.lg.km 
 -13.710551    1.279928 

              Estimate       S.E.    CI95.lo    CI95.hi Ztest      Zval
(Intercept) -13.710551 0.46745489 -14.626746 -12.794356   *** -29.33021
pop.lg.km     1.279928 0.05626785   1.169645   1.390211   ***  22.74705
Problem:
 Values of the covariate 'pop.lg.km' were NA or undefined at 0.57% (4 out of 
699) of the quadrature points
plot(effectfun(PPM1, "pop.lg.km", se.fit=TRUE), main=NULL, 
     las=1, legendargs=list(cex=0.8, xpd=TRUE, inset=c(1.01, 0) ))
Warning in min(D[scaledlegbox]): no non-missing arguments to min; returning Inf

Overlays

Overlays

  • Methods for identifying optimal site selection or suitability

  • Apply a common scale to diverse or disimilar outputs

Getting Started

  1. Define the problem.

  2. Break the problem into submodels.

  3. Determine significant layers.

  4. Reclassify or transform the data within a layer.

  5. Add or combine the layers.

  6. Verify

Boolean Overlays

  • Successive disqualification of areas

  • Series of “yes/no” questions

  • “Sieve” mapping

Boolean Overlays

  • Reclassifying

  • Which types of land are appropriate

nlcd <-  rast(system.file("raster/nlcd.tif", package = "spDataLarge"))
plot(nlcd)

Boolean Overlays

  • Which types of land are appropriate?
nlcd.segments <- segregate(nlcd)
names(nlcd.segments) <- levels(nlcd)[[1]][-1,2]
plot(nlcd.segments)

Boolean Overlays

  • Which types of land are appropriate?
srtm <- rast(system.file("raster/srtm.tif", package = "spDataLarge"))
slope <- terrain(srtm, v = "slope")
par(mfrow=c(1,2))
plot(srtm)
plot(slope)
par(mfrow=c(1,1))

Boolean Overlays

  • Make sure data is aligned!
suit.slope <- slope < 10
suit.landcov <- nlcd.segments["Shrubland"]
suit.slope.match <- project(suit.slope, suit.landcov)
suit <- suit.slope.match + suit.landcov

Boolean Overlays

par(mfrow=c(1,3))
plot(suit.slope)
plot(suit.landcov)
plot(suit)
par(mfrow=c(1,1))

Challenges with Boolean Overlays

  1. Assume relationships are really Boolean

  2. No measurement error

  3. Categorical measurements are known exactly

  4. Boundaries are well-represented

A more general approach

  • Define a favorability metric

\[ \begin{equation} F(\mathbf{s}) = \prod_{M=1}^{m}X_m(\mathbf{s}) \end{equation} \]

Weighted Linear Combinations

\[ \begin{equation} F(\mathbf{s}) = \frac{\sum_{M=1}^{m}w_mX_m}{\sum_{M=1}^{m}w_i} \end{equation} \]